#include <iostream>
#include<stdio.h>
#include<time.h>
#include<pthread.h>
#include <arm_neon.h>
using namespace std;
const int N=1000;
const int NUM_THREADS=7;
float **m ;
void serial(float **m) {
    for (int k = 0; k < N; k++)
    {
        for (int j = k + 1; j < N; j++)
            m[k][j] = m[k][j] / m[k][k];
        m[k][k] = 1.0;
        for (int i = k + 1; i < N; i++)
        {
            for (int j = k + 1; j < N; j++)
                m[i][j] = m[i][j] - m[i][k] * m[k][j];
            m[i][k] = 0;
        }
    }

}
void m_reset(float **m) {
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < i; j++)
            m[i][j] = 0;
        m[i][i] = 1.0;
        for (int j = i + 1; j < N; j++)
            m[i][j] = rand() % 100;
    }
    for (int k = 0; k < N; k++)
        for (int i = k + 1; i < N; i++)
            for (int j = 0; j < N; j++)
                m[i][j] += m[k][j];
}

typedef struct
{
    int k;//消去轮次
    int t_id;//线程id
}threadParam_t;

void *threadFunc_LU(void *param)
{
    threadParam_t *p=(threadParam_t*)param;
    int k=p->k;
    int t_id=p->t_id;
    int i=k+t_id+1;//获取自己的计算任务

    for(;i<N;i+=NUM_THREADS)
    {
        for(int j=k+1;j<N;j++){
            m[i][j] = m[i][j] - m[i][k] * m[k][j];
            }
        m[i][k] = 0;
    }
    pthread_exit(NULL);
    return NULL;
}
void *threadFunc_SIMD(void *param)
{
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k;
    int t_id = p->t_id;
    int i = k + t_id + 1;//获取自己的计算任务

    for (; i < N; i += NUM_THREADS)
    {
        float32x4_t vaik;
        vaik=vld1q_dup_f32(m[i]+k);
        int a=0;
        for (int j = k + 1; j +4< N; j+=4,a=j)
        {
            float32x4_t vakj=vld1q_f32(m[k]+j);
            float32x4_t vaij=vld1q_f32(m[i]+j);
            float32x4_t vx=vmulq_f32(vakj,vaik);
            vaij= vsubq_f32(vaij,vx);
            vst1q_f32(m[i]+j,vaij);
        }
        for(int j=a+1;j<N;j++)
        {
            m[i][j]=m[i][j]-m[k][j]*m[i][k];
        }
        m[i][k] = 0;
    }
    pthread_exit(NULL);
    return NULL;
    }
int main(){
    /********************普通单线程**********************/
    m = (float **)malloc(sizeof(float *) * N);
    for (int i = 0; i < N; i++) {
        m[i] = (float *)malloc(sizeof(float) * N);
}

    m_reset(m);
    struct timespec sts,ets;
    timespec_get(&sts,TIME_UTC);
    serial(m);
    timespec_get(&ets,TIME_UTC);
    time_t dsec=ets.tv_sec-sts.tv_sec;
    long dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0){
             dsec--;
    dnsec+=1000000000ll;
    }
    cout<<"LU分解单线程:";
    printf ("%lld.%09llds\n",dsec,dnsec);

    /********************LU分解多线程**********************/
    m_reset(m);
    timespec_get(&sts,TIME_UTC);
    for (int k = 0; k < N; k++)
    {//主线程做除法操作
        for (int j = k + 1; j < N; j++)
        {
            m[k][j] = m[k][j] / m[k][k];
        }
        m[k][k] = 1.0;
        //创建工作线程，进行消去操作
       
        pthread_t* handles = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));// 创建对应的 Handle
        threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));// 创建对应的线程数据结构
        //分配任务
        for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        {
            param[t_id].k = k;
            param[t_id].t_id = t_id;
        }
        //创建线程
        for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        {
            pthread_create(&handles[t_id], NULL, threadFunc_LU, (void*)&param[t_id]);
        }

        //主线程挂起等待所有的工作线程完成此轮消去工作
        for (int t_id = 0; t_id < NUM_THREADS; t_id++)
            pthread_join(handles[t_id], NULL);
        free(handles);
        free(param);
    }
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0){
             dsec--;
    dnsec+=1000000000ll;
    }
    cout<<"LU分解多线程:";
    printf ("%lld.%09llds\n",dsec,dnsec);

    /********************SIMD多线程**********************/
    m_reset(m);
    timespec_get(&sts,TIME_UTC);
     for(int k=0;k<N;k++)
    {
        float32x4_t vt=vld1q_dup_f32(m[k]+k);
        int a=0;
        for(int j=k+1;j+4<=N;j+=4,a=j)
        {
             float32x4_t va;
             va=vld1q_f32(m[k]+j);
             va=vdivq_f32(va,vt);
             vst1q_f32(m[k]+j,va);
        }
        for(int j=a+1;j<N;j++)
        {
            m[k][j]=m[k][j]/m[k][k];
        }
        m[k][k]=1.0;
    
        //创建工作线程，进行消去操作
        
        pthread_t* handles=(pthread_t*)malloc(NUM_THREADS*sizeof(pthread_t));// 创建对应的 Handle
        threadParam_t* param=(threadParam_t*)malloc(NUM_THREADS*sizeof(threadParam_t));// 创建对应的线程数据结构
        //分配任务
        for(int t_id=0;t_id<NUM_THREADS;t_id++){
            param[t_id].k=k;
            param[t_id].t_id=t_id;
        }
        //创建线程
        for(int t_id=0;t_id<NUM_THREADS;t_id++){
            pthread_create(&handles[t_id],NULL,threadFunc_SIMD,(void*)&param[t_id]);
        }

        //主线程挂起等待所有的工作线程完成此轮消去工作
        for(int t_id=0;t_id<NUM_THREADS;t_id++)
            pthread_join(handles[t_id],NULL);
        free(handles);
        free(param);
    }
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0){
             dsec--;
    dnsec+=1000000000ll;
    }
    cout<<"SIMD多线程:";
    printf ("%lld.%09llds\n",dsec,dnsec);
   
    

    return 0;
}